﻿#include "TroposphericCorr.h"
#include "literal.h"
#include <iomanip>

// 默认构造函数：使用标准大气参数
TroposphericCorr::TroposphericCorr()
{
    // 初始化标准大气参数
    stationTemperature = 288.15;    // 15°C转开尔文
    stationPressure = 1013.25;      // 标准大气压（毫巴）
    stationHumidity = 50.0;         // 50%相对湿度

    // 创建WGS84坐标系对象
    GPSConstant gpsConst;
    coordSys = new WGS84(gpsConst.get_refEllip());
}

// 带参数构造函数：设置气象参数
TroposphericCorr::TroposphericCorr(double temperature, double pressure,
    double humidity)
{
    // 设置气象参数
    stationTemperature = temperature;
    stationPressure = pressure;
    stationHumidity = humidity;

    // 创建WGS84坐标系对象
    GPSConstant gpsConst;
    coordSys = new WGS84(gpsConst.get_refEllip());
}

// 析构函数：释放动态分配的内存
TroposphericCorr::~TroposphericCorr()
{
    if (coordSys != nullptr)
    {
        delete coordSys;
        coordSys = nullptr;
    }
}

// 设置气象参数
void TroposphericCorr::set_meteorological_params(double temperature,
    double pressure,
    double humidity)
{
    stationTemperature = temperature;
    stationPressure = pressure;
    stationHumidity = humidity;
}

// Saastamoinen模型对流层延迟改正（重载1：使用坐标）
double TroposphericCorr::saastamoinen_corr(const XYZ& receiverPos,
    const XYZ& satPos)
{
    // 计算高度角
    double elevation = compute_elevation(receiverPos, satPos);

    // 获取接收机高程
    BLH receiverBLH = coordSys->XYZ2BLH(receiverPos);
    double height = receiverBLH.H;

    // 调用重载函数
    return saastamoinen_corr(elevation, height, stationTemperature,
        stationPressure, stationHumidity);
}

// Saastamoinen模型对流层延迟改正（重载2：使用参数）
double TroposphericCorr::saastamoinen_corr(double elevation, double height,
    double temperature, double pressure,
    double humidity)
{
    // 将高度角转换为弧度
    double elevRad = elevation * D2R;

    // 计算水蒸气压
    double waterVaporPressure = compute_water_vapor_pressure(temperature, humidity);

    // 计算干分量延迟
    double dryDelay = compute_dry_component(elevation, pressure, temperature, height);

    // 计算湿分量延迟
    double wetDelay = compute_wet_component(elevation, temperature, waterVaporPressure);

    // 总对流层延迟
    double totalDelay = dryDelay + wetDelay;

    return totalDelay;
}

// 计算高度角
double TroposphericCorr::compute_elevation(const XYZ& receiverPos,
    const XYZ& satPos)
{
    // 计算ENU坐标
    ENU enuPos = coordSys->XYZ2ENU(receiverPos, satPos);

    // 计算高度角（度）
    double range = sqrt(enuPos.e * enuPos.e + enuPos.n * enuPos.n + enuPos.u * enuPos.u);
    double elevation = asin(enuPos.u / range) / D2R;

    return elevation;
}

// 计算水蒸气压
double TroposphericCorr::compute_water_vapor_pressure(double temperature,
    double humidity)
{
    // 计算饱和水蒸气压（毫巴）
    double saturationPressure = 6.1078 * exp((17.269 * (temperature - 273.15)) /
        (temperature - 35.86));

    // 计算实际水蒸气压
    double waterVaporPressure = (humidity / 100.0) * saturationPressure;

    return waterVaporPressure;
}

// 计算干分量延迟
double TroposphericCorr::compute_dry_component(double elevation, double pressure,
    double temperature, double height)
{
    double elevRad = elevation * D2R;

    // 高度改正因子
    double heightCorr = 1.0 - (0.0026 * cos(2.0 * 0.0) - 0.00000028 * height);

    // 干分量天顶延迟
    double zenithDryDelay = (0.0022768 * pressure) /
        (1.0 - 0.00266 * cos(2.0 * 0.0) - 0.00000028 * height);

    // 干分量映射函数
    double dryMappingFunc = 1.0 / (sin(elevRad) + 0.00143 /
        (tan(elevRad) + 0.0445));

    // 干分量延迟
    double dryDelay = zenithDryDelay * dryMappingFunc;

    return dryDelay;
}

// 计算湿分量延迟
double TroposphericCorr::compute_wet_component(double elevation, double temperature,
    double waterVaporPressure)
{
    double elevRad = elevation * D2R;

    // 湿分量天顶延迟
    double zenithWetDelay = (0.002277 * (1255.0 / temperature + 0.05) *
        waterVaporPressure) / temperature;

    // 湿分量映射函数
    double wetMappingFunc = 1.0 / (sin(elevRad) + 0.00035 /
        (tan(elevRad) + 0.017));

    // 湿分量延迟
    double wetDelay = zenithWetDelay * wetMappingFunc;

    return wetDelay;
}

// 测试函数：验证对流层延迟改正功能
void TroposphericCorr::test(const XYZ& receiverPos, const XYZ& satPos)
{
    cout << "/* 对流层延迟误差改正测试 */" << endl;

    // 计算高度角
    double elevation = compute_elevation(receiverPos, satPos);

    cout << "卫星高度角: " << fixed << setprecision(2) << elevation << " 度" << endl;
    cout << "气象参数:" << endl;
    cout << "  温度: " << stationTemperature - 273.15 << " °C" << endl;
    cout << "  气压: " << stationPressure << " 毫巴" << endl;
    cout << "  湿度: " << stationHumidity << " %" << endl;

    // 计算对流层延迟改正
    double tropoCorr = saastamoinen_corr(receiverPos, satPos);

    cout << "Saastamoinen模型对流层延迟改正: " << scientific << setprecision(4)
        << tropoCorr << " 米" << endl;
    cout << endl;
}